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Abstract 

Marine ecosystem models are developed to understand and simulate the biogeo¬ 
chemical processes involved in marine ecosystems. Parekh, Follows and Boyle 
introduced the PO^-DOP-Fe model of the coupled phosphorus and iron cycles 
in 2005. Especially the part describing the phosphorus cycle (PO 4 -DOP model) 
is often applied in the context of parameter identification. The mathematical 
analysis presented in this study is concerned with the existence of solutions and 
the reconstruction of parameters from given data. Both are important questions 
in the numerical model’s assessment and validation not answered so far. In this 
study, we obtain transient, stationary and periodic solutions (steady annual 
cycles) of the PO^-DOP-Fe model equations after a slight change in the equa¬ 
tion modeling iron. This result confirms the validity of the solutions computed 
numerically. Furthermore, we present a calculation showing that four of the 
POi-DOP model’s parameters are possibly dependent, i.e. different parameter 
values might be associated with the same model output. Thereby, we identify a 
relevant source of uncertainty in parameter identification. On the basis of the 
results, possible ways to overcome this deficit can be proposed. In addition, the 
stated mathematical conditions for solvability are universal and thus applicable 
to the analysis of other ecosystem models as well. 
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identification 


1. Introduction 


Being a part of the global carbon cycle marine ecosystems considerably in¬ 
fluence on the earth’s climate. Marine ecosystem models, describing the bio¬ 
geochemical processes involved, provide an important tool to understand these 
processes and thus to predict the future concentration of carbon dioxide in 
oceans and atmosphere. 

One important exam ple for models of this kind is the PO^-DOP-Fe model 


Parekh et al. 


([200511 . describing the concentrations of three tracers phos- 


by 

phate (PO 4 ), disso lved organic phosphorus ( DOP ) and iron (Fe). Following 


Kriestetal 


(I2010I ), we use the synonymous name N-DOP-Fe model. The 


abbreviation N stands for the more general expression “nutrient”. 

Fixing the iron concentration, the evolving N-DOP model, consisting of 
the first two tracers, describes the marine phosphorus cycle. Despite t he low 


complexity and ass ociated simplifications, it is still in use (IKriest et al 


Parekh et al 


2012 


20061 1. An assessment on the basis of oceanic observations even 


indicates that, under certain c ircumstances, t he N-D OP model can compete 


with more complicated ones (Kriest et al 


2010 


2012). Due to the low com¬ 


plexity, it is often used for the purpose of testing numerical methods, e.g. in 


Priefietal, 


( 201311 . 


The third tracer expresses the influen ce o f iron on t he phosphorus cycle. Iron 


is an important nutrient. For instance, 


Kriest et al 


(12012 1 assume that their 


observed misfit of model output and observational data might be due to missing 
iron limitation. The model can also picture the changes in the phosphorus cycle 
induced by artificial fertilization with iron. 

The N-DOP-Fe model consists of three advection-diffusion-reaction equa¬ 
tions, each characterizing one tracer concentration on a three-dimensional ocean 
domain. The concentrations are influenced by ocean transport, i.e. advection 
and diffusion, and biogeochemical processes modeled by specific reaction terms. 
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The first part of this work studies the existence of transient, periodic and sta¬ 
tionary solutions of the N-DOP-Fe model equations. Periodic solutions (steady 
annual cycles) are the most relevant ones. Such results are desirable because the 
numerical model output used in applications is computed by approximating a 
solution of the original model equations. Therefore, the quality and reliability of 
the model output depends on the solvability of these equations. A comparable 
analysis of the N-DOP-Fe model equations has not been undertaken so far. 

Any model has to be calibrated, i.e. adapted to the ecosystem in ques¬ 
tion. This is achieved by the choice of the model’s parameters. Parameters, 
like e.g. remineralization rates or half saturation constants, are essential quan¬ 
tities characterizing the processes modeled. The N-DOP model includes seven 
parameters. 

Parameters can be estimated by means of laboratory experiments. If such 
measurements are too difficult or expensive, parameters are alternatively identi¬ 
fied via optimization with respect to observational data. Here, a parameter set 
minimizing the difference between model output and data is determined. An 
important area of application of the N-DOP model is the testing of nu meri¬ 


cal methods designed to solve such kind of minimization problems (jPriefi et al 


20131 ). The method in question is applied to synthetic data, i.e. data corre¬ 


sponding to known optimal parameters, assuming that a correct method is able 
to identify these. However, this assumption only holds true if all parameters 
are uniquely identifiable, i.e. if each possible model output is associated with a 
single parameter set. 

The second part of this work is dedicated to the question, unanswered so 
far, which of the N-DOP model parameters are uniquely identifiable. We name 
and justify different alternatives to alter the reaction terms in such a way that 
all parameters become identifiable. 

The mathematical results about the N-DOP-Fe model are introduced and 
explained to readers from other disciplines than mathematics. By outlining the 
mathematical proceeding, we intend to provide an overview about which math¬ 
ematical conditions and assumptions are responsible for the model’s properties. 
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These may be a guideline for the investigation and development of other models 
or alternative reaction terms. 

The paper is structured as follows: In the next section, we introduce the 
classical and the weak formulation of the N-DOP-Fe model equations, the lat¬ 
ter being the object of investigation. The next three sections each deal with one 
type of solution. The corresponding results concerning existence and, where 
possible, uniqueness are formulated and shortly justified. In Sec. [Gl we inves¬ 
tigate identifiability of parameters. In the last two sections, the results are 
discussed and conclusions are drawn. 


2 . Model equations 


The following introductio n of the N-D OP-Fe model is base d on the original 
work by 


Parekh et al. 

(2005 

) as well as on 

Roschat and Slawia 

(2014a 

t> 


The ecosystem is located in a three-dimensional bounded domain f2 deter¬ 
mined by the bounded water surface O' C R 2 and the depth h(x') > 0 at every 
surface point x' £ f 1'. The boundary T is the union of the surface T' and the 
boundary inside the water. 

The domain is separated into two layers, the euphotic, light-flooded zone fli 
up to a depth of h e := 120 m and the dark, aphotic zone beneath. The actual 
depth of the euphotic zone beneath some surface point x' is defined by h e {x') := 
min{h e , h(x')}. We split the surface into the part f V 2 := {x 1 £ £l']h{x') > h e } 
above the aphotic zone and the rest O'j := The boundary is analogously 

divided into the euphotic part Ti and the aphotic part T 2 . 

For the sake of a clearer terminology, we write henceforth 2/1 for the phos¬ 
phate concentration, y 2 for the concentration of DOP and 2/3 for the iron con¬ 
centration. The three tracers, assembled in the vector y = (2/1,2/2, IJ3), are 
characterized by the system of advection-diffusion-reaction equations 


d t yi + div(vyi) - div(KVi/i) - A y 2 + di(y) 
<9 t 2/2 + div(v 2 / 2 ) - div((cVj/ 2 ) + Aj/ 2 + d 2 ( y) 
<9(2/3 + div(v2/3) - div((cV 2 / 3 ) + JFe - A 2 / 2 i?Fe + d 3 (y) 


0 

0 ( 1 ) 
Sm, 
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each dependent on space and time in a finite time span [0, T], Regarding periodic 
solutions, a reasonable value for the final point of time T is one year. 

The second and third terms on the left-hand sides describe the ocean dynam¬ 
ics by means of the velocity field v and the diffusion coefficient n. Sometimes 
these te rms are sum marized in a linear operator or, if already discretized, in a 


matrix (Kriest et al. 


201C 


2012). To reduce computational effort, the models 


are often run in an “offline” mode, i.e. the influence the tracers have on the 
oceans dynamics is neglected. This has the advantage that the values of v and 
k can be precomputed with one run of an ocean circulation model. For that 
reason, v and k are assumed to be known in the theoretical analysis as well. 
Since turbulent exceeds molecular diffusion by far the values of k are assumed 
equal in all equations. 

Furthermore, since the ocean is a closed system, it is reasonable to assume 
that the velocity field v is divergence free and does not point out of the bound¬ 
aries. These two properties are crucial for the mathematical analysis. 

The reaction terms dj( y), Jj? e , Aj /2 and \y 2 RFe describe the biogeocliemical 
coupling. Therefore, they depend on one or several of the tracers. On the 
right-hand side of each equation, sources and sinks of the respective tracer are 
displayed. Only iron has a non-zero source term. 

In the following two subsections, we introduce and derive the reaction terms 
associated with the phosphorus cycle and the iron equation, respectively. Two 
further subsections deal with boundary conditions and weak solutions. 


2.1. The phosphorus cycle 

One important process is the remineralization of 2/2 into y \. In the first two 
equations of Eq. 0, it is modeled as a first order loss process with a remineral¬ 
ization rate A between 0 and 1. Being independent of light, this transformation 
takes place in the whole domain Q. The remaining processes, represented by 
the reaction terms d\ and c? 2 , differ according to the layers. In the euphotic 
zone, y\ is consumed by the photosynthesis of phytoplankton. This uptake is 
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modeled by 


G(yi,y 3 ) ■= a T 


V i 


2/3 


Jq—x^Kw 


\yx\ + K P \y 3 \ + K f \Ie *3K w \ + k i ' 

Here, the maximum uptake a > 0 is limited by the present concentrations of 
phosphate y\ and iron y 3 as well as insolation via Micliaelis-Menten kinetics. 
The corresponding saturation functions are equipped with the half saturation 
constants Kp, Kf, Kj. The absolute values in the denominators ensure that the 
fractions are mathematically well-defined because it is a priori unknown if the 
solutions 2/1 and y 3 are nonnegative. If this holds true, the above formulation 
is in accordance with Michaelis-Menten kinetics. Incidence of light is formal¬ 
ized by the insolation I = I(x',t) depending on the water surface and time. 
The incident light decreases exponentially with depth x 3 and the attenuation 
coeffici ent I\w f o r se awater. Below the euphotic zone, the values of / are zero 
(Paltridge and Platt, 1976 ). 

A fraction v € (0,1] of the uptake G is transformed into 2/2 while the rem¬ 
nants, integrated over the whole water column, are exported into the aphotic 
zone B 2 . The remineralization during the sinking of particles is described by a 
parameter b. 

The outlined processes are represented by nonlinear coupling terms with 
different appearance in each layer. In the euphotic zone Hi, they are given by 


di(y) := G(2/i,2/3) 
<k(y) : = ~vG(yi,y 3 ), 


and, in H 2 , by 

di(y) ■= -(1 - v) [ G(yi,y 3 )dx 3 ^- 

Jo h e \h e J 

<k(y) ■= 0 . 

We remark that the N-DOP model contains the seven parameters \,a,Kp, 
Ki , K w , 6 , v. 
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2.2. The iron equation 

In opposite to phosphorus, there is a source term for iron given by SVe := 
t 6Fi n which is non-zero only in the euphotic zone. The parameter /3 represents 
the solubility of iron in seawater and Fi n quantifies the aeolian source of iron. 

The reaction terms express how the phosphorus cycle, complexation and 
scavenging influence on the iron cycle. The iron concentration increases with re¬ 
mineralization and decreases with consumption of phosphorus. Being expressed 
in phosphorus units, these values are multiplied by the constant ratio I?p e in 
order to become iron units. Thus, iron increases with A yiF-Fz- The decrease 
induced by consumption is expressed by the coupling term 

d 3 (y) = G(y 1 ,y 3 )R Fe - 


2.2.1. Scavenging and complexation - original formulation 

The summand Jp e represents the influence of complexation and scavenging 
on the iron concentration. Since a certain amount of iron is complexed with 
organic ligand, the total iron y 3 is split into free iron Fe' and complexed iron 
FeL. Similarly, the total ligand Lt is the sum of free ligand L' and the com¬ 
plexed ligand, equal to the complexed iron FeL. These relatio ns are expre ssed 


by the formulae y 3 = Fe' + FeL and Lt = L' + FeL. 
Lt = 1 . 


Parekh et al 


(1200511 set 


Only the free iron is subject to scavenging which is thus modeled as the first 
order loss process Jp e := rk 0 CpFe'. Here, r, $ are numbers, k 3 is the initial 
scavenging rate and C p represents the particle concentration which decreases 
with depth. 

Since a feasible mathematical formulation requires that Jp e depends at 
least on onjy of the tracers yi,y 2 or y 3 , we express Fe ' using y 3 . To this 


end, 


Parekh et al 


({20051 1 additionally provide the equilibrium relationship K = 


FeL/Fe'L' with a positive constant K. Inserting the equivalent expression 
L' = FeL/KFe' into the equation for ligand we obtain 


Lt = FeL -f- L' = FeL - 


FeL 

KFe' 


= FeL 1 + 


KFe' 
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FeL can be replaced by y 3 — Fe'. This gives 


L T = (y 3 -Fe') (1 + ^ 7 ) 


= 2/3 ~Fe' + 


2/3 


1 

K' 


KFe' 

With the abbreviation H(y 3 ) := Lt + 1/A' — y 3 , this proves equivalent to 

Fe' 2 + H(y 3 )Fe'-^ = 0. 

The two solutions of this quadratic equation are known. It can be easily shown 
that one of them has only negative values and is therefore inappropriate to 
describe the amount of free iron. Thus, we find 


Fe\y 3 ) = ~H(y 3 ) + 


H(yzf 


2/3 

K' 


To ensure that the square root is real we show that the radicand r := (Lt + 
1/K — 2/ 3 ) 2 /4 + 2 / 3 /A' is positive. This is obvious whenever y 3 is positive, since 
in this case r > 2 / 3 /A' > 0. In the non-positive case, we transform r into 

2 


Lt + — ~ 2/3 


2/3 

K 


= - I Lt + 


— - ( Lt + 


I< 


2/3 


r?/! 


2/3 

K 


1 


= 7 ( Lt + % ) + 7 % - o 2/3 ( L t - — 


1 


K 


Since K tends to be a large number and Lt takes a value around one (jParekh et al 


( 20051 ) set A' = exp(ll) and Lt = 1) the expression Lt — 1/K is nonnega¬ 
tive. Therefore, y 3 < 0 yields —y 3 (LT — l/K)/2 > 0 and, as a consequence, 
r > (Lt + l/A') 2 /4 > 0. In particular, this result implicates that Fe' is differ¬ 
entiable everywhere on R. It can be shown that the derivative is positive and 
thus that Fe' is a monotonically increasing real function. 

In Fig. m we see that, as long as approximately y 3 < Lt = 1, the values for 
Fe' increase very slowly, for y\ > Lt they increase with a gradie nt of almost 


one. This behavior corresponds to the statement of 


Parekh et al. 


( 20051) that 


they “rapidly precipitate Fe' when FeT > At” meaning that, as long as ligand 
is available, only a small amount of iron remains free. 
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Figure 1: Graph of free iron Fe' in dependence of total available iron 1/3 assuming K = exp(ll) 


and Lt = 1 I Parekh et al 


200511 . 


Thus, the reaction term associated with scavenging and complexation de¬ 
pending on the third tracer 2/3 is determined by 

Jf e( 2 / 3 ) := rk 0 CpFe'(y 3 ). 

2.2.2. Scavenging and complexation - adjusted formulation 

The reaction term J Fein’s) derived above lacks one important property with 
regard to solvability. This is due to the fact that the function 

Fe ’ M = 4 [ Lt+ i ~ ■») + \j 1 ( Lt+ i -») +1 

modeling free iron levels off fast in negative direction. For that reason, we 
propose an alternative term FeF{y 3 ) for free iron. Fig. |T| indicates that Fe 1 
resembles a straight line with a slope of 1 for approximately y 3 > Lt- If y 3 < Lt, 
the curve tends to zero very fast. As a substitute for Fe' we therefore define a 
piecewise linear function, composed of a line with a slope of 1 and another line 
through zero. As their “meeting point” we determine ( LT,Fe'{LT )) where the 
function Fe' apparently turns upwards. Formally, the function FeF depending 
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Figure 2: Graph of free iron Fe' (red line) in comparison to t he alternative funct ion FeF 
(blue line) around ys = Lt assuming K = exp(ll) and Lt = 1 (jParekh et al.l . l2005h . 


on total iron 7/3 is defined by 


FeF(y 3 ) = 


2/3 + Fe'(L T ) 

^ 2/3 


Lt if 2/3 > Lt, 

if 2/3 < Lt- 


Figure [2] shows the difference between exemplary curves of Fe' and FeF 
around 2/3 = Lt- Comparing the complete curves reveals that, in the positive 
region, the piecewise linear function FeF lies slightly above the original and their 
distance remains small. In the negative region, the curves behave conversely. 
The distance increases with decreasing values of 2/3 because FeF has a constant 
slope and Fe' a decreasing one. However, negative values for 2/3 are not relevant 
in most applications. 

Taking into account these considerations, a possible alternative for Jp e might 
be the adjusted reaction term 

J{ 2/3) -=TkoCpFeF(y 3 ). 
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2.3. Boundary conditions 

The original N-DOP-Fe model formulation lacks explicit statements about 
the tracers’ behavior on the boundary. As this information is needed with 
respect to a weak formulation of the model equations (cf. Sec. 12.41 . suitable 
conditions are derived in this section. In general, a boundary condition of 
Neumann type for the j-th tracer has the form 


V,i )j ■ (Kr]) + bj( y) = 0 . 


( 2 ) 


The symbol 77 stands for the outward-pointing unit normal vector field on the 
boundary. Therefore, • (k rj) can be understood as the change of concentra¬ 

tion alongside the “conormal” vector nr). The term bj( y) describes the tracer 
coupling on the boundary. 

As to boundar y con ditions b ±, b 2 for the two equations of the phosphorus cy- 


Parekhetal 


(12005 1 indicate that “any remaining particulate organic matter 


cle, 

that reaches the bottom of the model domain is instantly remineralized”. Ad¬ 
ditionally, the model equations contain no sources or sinks for phosphate and 
DOP. This means that the total concentration (or mass) of the first two tracers 
remains constant with respect to time. We formalize the time-dependent total 
mass by 

mass(yi,y 2 ) := / (yi + V 2 )dx. (3) 

J n 

Mathematically, conservation of mass thus corresponds to the condition 

(1 f 

— mass(yi,y 2 ) = / dt{yi + U 2 )dx = 0 for all t. 

Jq 

On this basis we are able to derive boundary conditions. With the help of the 
model equations we replace the sum of the partial derivatives with respect to 
time dt(y 1 + y 2 ) and obtain 


0 = / div(v(yi + y 2 ) - n\7(y 1 + y 2 ))dx 
Jn 


(4) 


-6-1 


+ [ (1 - v)G(y 1 ,y 3 )dx - f (1 - v) [ Giyx.y^dx^^- f ^ ) dx. 

Jih Jo.2 Jo \ n e 
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Referring to the hindmost term as M, we obtain by inserting the definition of 
fl 2 and solving the emerging integral with respect to x 3 analytically 


M : = (1 - v) [ [ he G( yi ,y 3 )dx 3 t Yf^) & \x 3 dx' 

Jo'J0 "e Jli, V'W 

= G ^ Vl,y3 ^ dX3 {^~) [ _:C 3 6 ] } ^ 

= (1 - V) J^G( yi ,y 3 )dx 3 ^1 - j dx'. 

Taking into account the definition of fii, the integral with respect to Hi in 
Eq. (Q} equals 


/ (1 - v)G( yil y 3 )dx = / / (1 - v)G{y i,y 3 )dx 3 dx'. 

JJJ 0 

Finally, Gaufi’ divergence theorem yields for the first integral in Eq. ® 

/ div(v(yi + y 2 ) - nS7(yi + y 2 ))dx 
J n 

= / (v • 77(2/1 + 2/2) - V(2/1 + 2/2) • (m]))ds = - / V(2/1 + 2/2) • (K 7 j)ds. 


J r ./r 

The last equality sign holds because of the assumptions about v. Combining 
the results obtained and heeding SI' = UH^, we transform Eq. into 


V(yi +y 2 ) • (KT])ds 


= (1 - ^) 


/ni Jo 


G{y\,y 3 )dx 3 dx' + 


nt Jo 


G(yi,y 3 )dx 3 



Since the integral over S7' corresponds to the boundary integral over T^, the 
last result allows us to formulate proper boundary conditions. First of all, the 
right-hand side is independent of y 2 such that a reasonable definition is b 2 (y) = 
0 on T. In accordance with the integrands of the above integral condition, we 
define the boundary coupling term for the first equation by 


My) : = 


-(! - v) Jo‘ G(yi,y 3 )dx 3 


in Ti, 


( 1 -^)Jo e G{yi,y 3 )dx 3 (^) in E 


as well as b\ (y) = 0 on T'. 
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The third tracer’s behavior on the boundary remains unspecified bv IParekh et al 


( 2005T) . Since there is a source of iron, it is not appropriate to claim conservation 
of mass. For this reason, the boundary conditions are limited to expressing that 
iron does not escape through the boundary. Speaking in mathematical terms, 
we impose the homogeneous boundary condition 63 (y) = 0 on T. 

2-4■ A weak formulation 

The analysis in this study refers to weak solutions. These solve a “weak 
formulation”, derived from the classical model equations (cf. Eq. (fljll by re¬ 
laxing the regularity requirements. Weak formulations are often preferred to 
classical ones because they are defined on Hilbert spaces and thus a well in¬ 
vestigated theory about existence and uniqueness of solutions is available. The 
consideration of weak solutions is justified by the fact that every weak solution 
automatically solves the classical formulation, provided that the corresponding 
regularity assumptions are fulfilled. 

An important area of application is optimal control theory which is designed 
to solve optimality problems governed by weak formulations of differential equa¬ 
tions. Since parameter identification leads to such a problem, weak solutions 
will be required in the associated Sec. [G] 

A weak formulation is derived by integrating the original equation multiplied 
with a “test function”. By integration by parts, regularity is partly transferred 
to the test function. An appropriate weak formulation for the N-DOP-Fe 
model equations is 

[ {{y[,wi) + B(y 1 ,vui)+ [ (~\y 2 + di(y))widx + f b 1 (y)w 1 ds}dt = 0 

Jo Jn Jr 

[ {( 2 / 2 ^ 2 ) + B(y 2 ,w 2 ) + [ (Xy 2 + d 2 (y))w 2 dx}dt = 0 

Jo J n 

/ {( 2 / 31 ^ 3 ) + B(y 3 ,w 3 ) + / (JFe(y 3 ) - Xy 2 R Fe + d 3 (y))w 3 dx}dt 


= / SFeWodxdt. 

Jo Jn 

By wi,w 2 ,W 3 we denote the test functions. The angle brackets indicate that 
the temporal derivatives j/' are elements of a “dual space”, i.e. maps on the 
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space of test fu nctions. This kind of differentiability, usually referred to as 
“distributional” (IZeidlerl . Il990l Def. 23.15), is weaker than differentiability in 
the usual sense. The time-dependent term B , defined by 

B(yj,Wj) := / (nVi/j ■ S7vjj)dx + / di v(vyj)wjdx 
Jq J q 

summarizes the ocean transport. The first summand is the resulh_of applying 
integration by parts or, more precisely, Green’s identity (IStraussl . 120081 Chap. 7, 
Eq. (Gl)) to the diffusion term. In this connection, the boundary conditions, 
formerly formulated on their own, enter the weak formulation. 

Leaving out test functions and integrals, the weak formulation becomes the 
equivalent “operator equation” 


y'l + B{yi) - \y 2 = -di(y) - 6i(y) 
2/2 + B (y 2 ) + Xy 2 = -d 2 (y) 

2/3 + B (y 3 ) + Jf e( 2 / 3 ) ^ A y 2 R Fe = S Fe - d 3 ( y), 


(5) 


valid in the already mentioned dual space. In the following, we deal with the 
adjusted version of Eq. © distinguished by the appearance of J instead of J Fe 
in the last equation. 

The weak formulation introduced here is derived in detail bv lRoschat and Slav 
( 2014al) . Further inf ormation about weak formulation s and o perator equati ons 
can be found e.g. in lTroltzsch (2010|); [Gajewski et al . (1197411 : 


Zeidle r (19901 ). 


Three different types of solution (transient, periodic, stationary) will be dealt 
with successively in the following sections. 


3 . Transient solutions 

This section deals with solutions with a given initial concentration. Be¬ 
side the weak formulation, they satisfy the initial value condition yj(x, 0) = 
Uq(x) for x £ 12 with a prescribed concentration y° 0 at the point of time t = 0. 
Solutions of this kind are called transient. 
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The adjusted N-DOP-Fe model equations have a unique transient solution 
for every at least quadratically integrable initial concentration. Furthermore, 
the solution depends continuously on the initial concentration. 

To justify this assertion we refer to our recent analysis of system s a nalog ous 
to Eq. © with an arbitrary number of equations (IRoschat and Slawig, 2014at ). 


The main challenge is the t reatm ent of the nonlinear coupling terms. In our 
analysis, we followed [Evan s (119981 Sect. 9.2, Thm. 2) by converting the problem 
of solving the weak formulation into a fixed point problem. To achieve this, an 
arbitrary fixed vector z = (zi, Z 2 , z 3 ) is inserted into the reaction terms leading 
to the problem 


y[ + B{y 1 ) = Xz 2 - di(z) - bi(z) 
y'2 + B(y 2 ) = -Xz 2 - d 2 (z) 

2/3 + B{y 3 ) + J{y 3 ) = $Fe + Xz 2 R,Fe - d 3 (z) 

Vj{ 0) = Vo for all j <E {1, 2, 3}. 

Unique solvab i lity of linear or mono t one equations is a well -investigated problem 


(Zeidler, 

1990: 

Ladvzenskava et al.. 

1968; 

Troltzsch. 

2010) 


2010(1 . While the first two 


equations are linear, the third is monotone because the reaction term J is based 
on a monotonically increasing real function (cf. Sec. 12.211 . In addition, both 
B and J have a linear growth and B is monotone and linear. Therefore, the 
system depending on z has a unique solution y. 

Obviously, a fixed point of the map z 4 y = (y i, y 2 ,y 3 ) corre sponds to a 


transient solution. Banach’s Fixed Point Theorem (Zeidler 


19861 Thm. l.A) 


yields the existence of a unique fixed point provided that the reaction terms 
on the right-hand side are Lipschitz continuous with a Lipschitz constant in¬ 
dependent of time. In the N-DOP-Fe model, these terms are either linear 
(multiplication with A) or based on the uptake function G. However, G proves 
Lipschitz continuous on R 2 , equipped with an arbitrary norm ||.||, according to 
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the estimation 


< a 


G(wi,u) 3 ) < a 

2 /i 

2/3 

w 1 ui 3 

I2/1I + Kp I2/3I + K F 

|uq + Kp IC 3 + Kp 

( 

2 /i 

Wi 

2/3 



\yi\+K P wi| + Kp 

I2/3 1 + K F 



w 1 


|wi| + K P 


2/3 


W 3 


|y 3 | + K F \w 3 \ + K F 


< a 


\y i ~ Wl \ + \V 3 ~w 3 \ ) < C\\(yi,y 3 ) - (wi,w 3 )||, 


using that saturation functions are Lipschitz continuous and that the Lipschitz 
constant equals the inverse of the half saturation constant. Furthermore, the 
absolute value of the saturation functions is bounded by one. C is a constant 
depending on a, Kp, Kp. 

It can be stated as a general observation that model equations have transient 
solutions if all reaction terms are either Lipschitz continuous with a Lipschitz 
constant independent of time or if they are monotone with linear growth. 


4 . Periodic solutions 

In this section, we consider the existence of steady annual cycles or peri¬ 
odic solutions. These are characterized by the additional condition yj(x, 0) = 
yj(x,T) for x G 17 signifying that the solution reaches its initial value again at 
the time T. Interpreting T as one year, periodic solutions correspond to steady 
annual cycles, generally approximated via spin-up or fixed-point iteration. Un¬ 
fortunately, it is seldom possible to prove analytically, that the fixed-point iter¬ 
ation converges. 

Given a fixed C > 0, there is a periodic solution of the adjusted N-DOP- 
Fe model with the additional property that the sum of the first two tracers’ 
mass (hereafter called total mass although iron is excluded) is equal to C. In 
particular, a periodic solution is not unique. The time-dependent total mass is 
formalized in Eq. ©. Remark that, due to the choice of reaction terms and 
boundary conditions, every solution has a constant total mass with respect to 
time. 
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The justification of the existence result bases on our previous study about 


periodi c so lvability of the N-D OP model, assuming a constant distribution of 


iron (Roschat and Slawid . 


2014bj). We will present our proceeding and indicate 


how to extend the argumentation to the model with three equations. 

In opposite to the transient case, the model equations are simplified by 
inserting a fixed z only into the reaction terms dj and bj. In a first step, we 
solve the periodic problem 


Vi + B (yi) ~ Aj/2 = ~di(z) - b±(z) 

2/2 + B ( yz ) + ^2/2 = ~ d 2 ( z ) 

2/3 + b (V 3) + J{y3 ) - A y 2 R Fe = S Fe - d 3 (z) 

Vj (0) = Vj(T) for all j G{ 1,2,3} 


and find a fixed point of the map zH-y afterwards. Since there is no standard 
technique for periodic solvability of such systems, we had to develop a method 
of our own. Its presentation provides insight into the behavior of closed systems. 

Being partly decoupled, the three equations of the system above can be 
solved consecutively. The second equation has a unique periodic solution de¬ 


pending on z because the linear summand +Xy 2 allows the estimate 

/ [ Xy 2 y 2 dxdt > X\\y 2 \\ 2 


( 6 ) 


where ||.|| is the norm in a Hilbert space of functions depending on space and 
time. This condition causes the corresponding fixed-point iteration to converge. 
Since the first equation lacks a comparable summand, we first solve S' + B(S) = 
— di(z) — 61 (z) — d 2 (z), the equation for the sum S = yi + y 2 . A standard 
existence theorem provides a uniq ue so lution in a s pecia l function space with the 
additional condition J n Sdx = 0 ( Gaiewski et all 11974 Chap. VI, Thm. 1.4). 
The properly chosen boundary conditions (and therefore conservation of mass) 
ensure that S belongs to the desired function space. Due to the choice of S, the 
function y\ = S — y 2 + |ST|solves the first equation. Here, the symbol |fi| 
stands for the measure of the domain SI. The solution y = (yi,y 2 ) is unique 
and satisfies mass(yi, y 2 ) = f n \Vl\~ l Cdx = C. 
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The third equation can be solved using the mentioned standard theorem as 
well. The adjusted coupling term J is strictly monotone because FeF has a 
strictly increasing slope (cf. Sec. 12 . 21 ) and all other factors are positive. Addi¬ 
tionally assuming a positive lower threshold c p > 0 for the particle concentration 


C p , the term J satisfies 


(7) 


/ / J{yi)y^dxdt > C1II2/3II 2 — c *21 |2/3II- 

Jo Jn 

This condition is slightly weaker than the one displayed in Eq. © but still 
sufficient to ensure that the operator B + J is “coercive” which is a crucial 
property in the existence theorem. Actually, the original reaction term Jp e was 
replaced because it lacks coercivity. Taking into account the properties of J, 
the standard theorem yields a unique periodic solution y 3 depending on z. 

As in the transient case, a fixed point of the map z 1 —> y = (2/1,2/2 5 2/3) 
corresponds to a pe riod ic sol ution of the original problem. Schauder’s Fixed 


Point Theorem (Zeidler 


19861 Thm. 2.A) yields the desired fixed point because 


the coupling terms dj and bj are bounded independently of z. This is due to 
the uptake G(z±, z 3 ) whose absolute value is bounded by 1 . 

As a general observation, we can state that a proof of periodic solvability 
requires a specialized proceeding adapted to the model in question. Concerning 
the choice of reaction terms, it can be said that every function bounded inde¬ 
pendently of the inserted argument is allowed. Furthermore, the j-th equation 
can (and to a certain extent has to) contain reaction terms which are monotone 
and coercive, i.e. satisfy an estimate analogous to Eq. 0 , and depend only on 
the j -th tracer. 


5. Stationary solutions 

Stationary solutions represent the constant tracer concentrations reached 
with a fixed forcing. Consequently, they solve a time-independent variant of 
the original model Eq. 0 , characterized by vanishing temporal derivatives and 
temporally constant terms v, n, dj ,bj, Jp e - The corresponding adjusted weak 
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formulation including test functions is 


B{yi,w 1 )+ / (-Xy 2 +d^y^wxdx + / b 1 {y)w 1 ds = 0 
Jn Jr 

B(y 2 ,w 2 ) + / {Xy 2 + d 2 (y))w 2 dx = 0 (8) 

Jn 

B(y 3 ,w 3 ) + / {J(y 3 ) - Xy 2 R Fe + d 3 (y))w 3 dx = / S Fe w 3 dx. 

Jn Jn 

Since the equation is independent of time, the integrals over [0,T] are missing. 
For the same reason, the test functions, like the solutions, only depend on the 
spatial coordinates. 

It is important to bear in mind that stationary solutions are not periodic 
in the sense of the last section although they are constant with respect to time 
and thus initial and terminal values coincide. The reason is that stationary 
solutions in the sense of this study correspond to a constant forcing, while pe¬ 
riodic solutions solve the equations with temporally variable summands. Thus, 
the periodic solutions of the last section are not necessarily constant. 

The adjusted N-DOP-Fe model equations have a stationary solution, i.e. a 
solution of Eq. ©, for every total mass C > 0. 

An analysis of stationary solutions has not been published yet. However, 
since they are closely related to periodic solutions, the considerations of the 
last section can be transferred. As in the periodic case, the system in Eq. m 
is simplified by inserting a fixed z into the reaction terms dj and bj. The 
equations are solved in the same orde r as the period ic ones, now with the help 
of the theorem of Browder and Minty flZeidlen . Il989l Thm. 26.A) instead of the 
periodic standard theorem. The theorems are closely related. Both claim the 
assumptions continuity, monotonicity and coercivity all of which have already 
been investigated in the last section. Finally, Schauder’s Fixed Point Theorem 
is applied in exactly the same way as before. 

Due to the analogous proceeding, the conclusions concerning the choice of 
reaction terms are the same as in the last section. 
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6 . Parameters in the N-DOP model 


Parameter identification is one of the most important areas of application 
of the N-DOP model. As explained in the introduction, parameters are iden¬ 
tified via minimizing the difference between model output, depending on the 
parameters, and original data. Numerical optimization methods to be tested 
are applied to synthetic data and judged by their ability to identify the corre¬ 
sponding known optimal parameters. 

A criterion for the explanatory power of this approach is the parameters’ 
unique identifiability. Parameters are called uniquely identifiable if each set 
of parameter values is associated with a single model output. Otherwise, the 
parameters are called dependent. 

Thus, information about dependencies in the N-DOP model contributes to 
the validation of numerical tests and the interpretation of results like those of 
Priefi et all (120131! whose method did not identify all parameters correctly. 


6.1. Investigation of identifiable and dependent parameters 

In this section, we investigate the seven N-DOP model parameters as to 
possible dependencies. As a result, we will be able to tell which parameters 
are uniquely identifiable and thus suited for testing purposes and which are 
(supposedly) dependent. 

To this end, we assume that the equations associated with the parameter sets 
ui = (\i,ai,Kpi,Kii,Kwi,bi,vi) and u 2 = (A 2 , a 2 , Kp 2 , Kj 2 , Kw 2 , b- 2 , ^ 2 ) 
have the same nontrivial solution y = (yi,y 2 ). It is not relevant if the solution 
is periodic, transient or stationary. Parameters that prove equal in both param¬ 
eter vectors are uniquely identifiable. We limit the investigation to the natural 
case that y has two nontrivial components because otherwise there are obvious 
dependencies. For instance, the solution y\ = 0 allows arbitrary parameters 
a, Kp, Kj, Kyy and y 2 = 0 allows an arbitrary A. 

It is reasonable to restrict the considerations to positive a, v < 1 and K[ > 0. 
The value a = 0 is not a likely maximum production and obviously effects 
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dependencies since a = 0 and Kj. Kp. Kw arbitrarily chosen lead to the same 
solution. In case v = 1, the export into the aphotic zone is zero because all 
consumed phosphate is transformed into DOP and thus the sinking parameter 
b can be chosen arbitrarily. Finally, a vanishing half saturation constant Kj = 
0 eliminates the influence of light. All other parameters are assumed to be 
nonnegative. 

The choice of y signifies that 

y[ + B (yi) = \yi - di(ui,y) - M u *,y) 

2/2 + B (y2) = -\y2 - d 2 (u tl y) 

hold for both i = 1,2. This time, we explicitly indicate the reaction terms’ 
dependence on the parameter vectors. Clearly, the left-hand sides are equal for 
both i and therefore also the right-hand sides, i.e. 

A 12/2 - di(ui,y) - 6i(ui,y) = X 2 y 2 ~ di(u 2 ,y) - &i(u 2 ,y) (9) 

-Ai 2 / 2 - d 2 (ui,y) = -A 2 y 2 - d 2 (u 2 ,y). (10) 

In order to draw conclusions about the parameters, we utilize the reaction terms’ 
specifications in each part of (see Sec. 12.11) . Equation (fTUl) . considered in fl 2 , 
reveals that Ai y 2 = \ 2 y 2 . Since y 2 is not trivial, this shows Ai = A 2 , i.e. unique 
identifiability of A. 

Taking into account Ai = A 2 , Eq. © in fli yields 

G(ui,yi) = G(u 2 ,j/i). (11) 

From Eq. (1113) in fli, we conclude with the help of Eq. (HD 

0 = zqG(ui,yi) - i/ 2 G(u 2 ,yi) = {ui - v 2 )G( ui,yi) + v 2 (G(vli, yi) - G(u 2 ,yi)) 

= {v-l - v 2 )G( ui,yi). 

Since a > 0 and j/i is not trivial by assumption we obtain G(u 1 ,y 1 ) ^ 0 and 
thus vi = v 2 , i.e. unique identifiability of v. 

In the next step, we deal with the parameter b. Since otherwise this param¬ 
eter would not occur in the model equations, we assume that the integral over 
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G(ui, yi) with respect to depth is not zero everywhere. Taking into account the 
aphotic boundary T 2 , it follows from Eq. © 


rh e 

0= / G(ui, 2 /i)dx 3 
Jo 



G(ui,yi)dx3 


(my b '- £ G{ ^ l)dx , 




—b 2 


+ (yjp-) b2 [\ G ^yi)-G{n 2 ,y r))dx 3 

due to 1 — v ^ 0. The last summand vanishes according to Eq. m Because 
of the assumption concerning the integral over G(ui,yi), we conclude 


feQrOy* 1 fh(x') y b2 Q _ 

he ) \ h e J 

The fraction q := h(x')/h e is strictly less than 1 for at least one x' since otherwise 
the domain would lack the aphotic part and the corresponding boundary. Since 
the natural logarithm In : R>o —>• R is bijective we obtain 


-hHQ) = Hq~ bl ) = Hq~ b2 ) = -& 2 Hq)- 


As ensured before, ln(g) / 0. It follows b\ = 62 . 

The remaining parameters are incorporated in the uptake function G. Equa¬ 
tion m states that 

yi j g -x 3 k w1 yi Ie -x 3 K W2 

Q 1 -— ----—-—-— Qjq ■ ■ - ~ 0. 

\yi \ + Kpi \Ie~ x 3 K wi \ + Kn |j/i | + Kp 2 \Ie~ X3Kw2 \ + K I2 

This equality provides the possibility to derive a condition for dependencies. 

The computation provided in App. |Appcndix A| shows that dependencies exist 

if and only if it is possible to depict the solution of the first equation by 

c\I + c 2 e X3C7 — C 3 e X3Cs 
12/11 _ al - c 5 e X3 ° 7 + c 6 e X3Cs 

with constants C\,... ,Cg satisfying the additional condition c\ = C2/C5 — (04 + 
l)c 3 /eg. In other words, if the absolute value of y± can be expressed in the 
indicated manner, it is possible to find two parameter vectors Ui,U 2 leading 
to the same solution. We cannot rule out that this condition holds for some 
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coefficients ci,..., eg. Thus, we have to assume that dependencies exist and the 
four parameters a, Kp, Kj, Kw may not be uni quel y id entifia ble. 


This result coincides with the finding of 


PrieBetal, 


(]2013f ) who were able to 


identify A, b and is (which they call a) very well. The results for the four other 
parameters were less satisfying. Whereas the optimal a was almost reached, the 
approximation hardly moved towards the optimal Kj. Concerning the other 
parameters, the accuracy of approximation lay in between. 


6.2. Elimination of dependencies 

To be valid and significant, tests of numerical methods should rely on models 
containing only identifiable parameters. In this section, we therefore propose 
two possible ways to eliminate the N-DOP model’s dependencies. 

The first possibility is to identify a reduced number of parameters after 
fixing the remaining ones. The fixed values have to be quantified in another 
way (experiments, estimations). In case of the N-DOP model we propose to 
fix A'/ and Kw- Then the five parameters A, a, Kp,b, is remain to be identified 
via optimization. 

We justify this suggestion by proving the unique identifiability of the five 
variable parameters. This property has already been shown for A, 6, is. Con¬ 
cerning the two remaining parameters, the analog of Eq. m 

2/i 2/1 

ai \yi\ + K P1 ~ a2 \ yi \ + K P2 ' 
can be simply transformed into 

oq _ \yi \ + Kpi _ i Kpi — Kp 2 

ot 2 \yi\ + Kp 2 \yi\+Kp 2 

The left-hand side is clearly constant. Since the only constant solution yi = 0 is 
excluded, the right-hand side is constant if and only if the last fraction vanishes, 
i.e. if Kpi = Kp 2 . We directly conclude ai/a 2 = 1, he. oq = a 2 . Thus, both 
remaining parameters are identifiable. 

We state that slightly reducing the number of unknown parameters in the 
N-DOP model eliminates all dependencies. Thus, the N-DOP model with the 
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Figure 3: Comparison of saturation function (blue) with a = K = 1 and the arc tangent 
scaled with the parameter /3 = 0.6 (red). 

reduced parameter set is recommendable to test and validate numerical methods 
reliably. 

A second possibility to construct a model with identifiable parameters is to 
replace the reaction terms being responsible for dependencies. This may involve 
a reduced number of parameters which are not immediately interpretable with 
respect to the biochemical processes modeled. 

In case of the N-DOP model, the results indicate that dependencies are 
caused by the product of saturation functions. A possible alternative for the 
saturation function ayi/(\yi\ + K) with two parameters a, K could be the scaled 
arc tangent /3 arctan(yi) with only one parameter /?. As we see in Fig. [3] the 
curves of both functions with appropriate parameter values behave similarly. 

Alternatives to the arc tangent are for example an approximation via the 
Taylor series or a linear function. Which one is suited best, depends e.g. on the 
values the solution y\ reaches. 
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7 . Discussion 


A numerical treatment of partial differential equations is usually preceded by 
a mathematical analysis. Numerical methods are designed to approximate exact 
solutions. Thus, if there is no exact solution, the numerically obtained results 
are dubious and unpredictable. A mathematical analysis therefore contributes 
to the validation and assessment of ecosystem models. In thi s work , we provide 


Parekh et al. (2005) 


correspondent information about the N-DOP model by 
which is extensively used. 

In the first part, the model formulation is stated in full mathematical detail. 
We additionally develop boundary conditions assuming conservation of mass. 
This condition is both common regarding ecosystem models and crucial in the 
proofs of periodic and stationary solvability. Furthermore, we explicitly derive 
the reaction term modeling complexation and scavenging of iron and find an 
alternative formulation with a similar behavior and convenient mathematical 
properties. Finally, a weak formulation of the model equations is specified. The 
exact mathematical formulation developed here is an essential premise for the 
following analysis. 

As a result, the analysis yielded that the adjusted N-DOP-Fe model and 
therefore, in particular, the (original) N-DOP model has transient, stationary 
and periodic solutions. 

Periodic solutions or steady annual cycles, characterized by equal initial and 
terminal values (e.g. at the beginning and the end of one year) are of particular 
interest as they are required in most applications. The analysis of periodic so¬ 
lutions provides an interesting insight into systems conserving mass. Assuming 
conservation of mass, which is commonly done in the context of ecosystem mod¬ 
els, we were able to show that there exists a distinct periodic solution for each 
mass contained in the phosphorus cycle. In particular, there are other periodic 
solutions than the trivial y = (0,0,0), corresponding to the mass zero. This 
result confirms the observation made in spin-up computations that every initial 
mass leads to a corresponding periodic solution. 
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As a conclusion, the results suggest that it is possible and meaningful to 
solve the N-DOP model equations numerically. 

Furthermore, the mathematical conditions found during the analysis enables 
the identification and elimination of problematic reaction terms. The reaction 
term for scavenging and complexation of iron for instance, could be sub¬ 
stituted by a suitable alternative J. 

In the second part, we investigated the N-DOP model’s parameters. Four of 
the seven parameters proved probably dependent and thus not uniquely identi¬ 
fiable. This result suggests that unsatisfying results in parameter identification 
might be due to the model itself instead of an inadequate numerical method. As 
a consequence, the original N-DOP model seems unsuited for tests of parameter 
identification methods. 

The analysis reveals that dependencies originate from the modeling of the 
biological uptake limited by iron, phosphate and light. Thus, a model with¬ 
out dependencies either lacks this term altogether or has a reduced number of 
parameters to be identified. We showed that one possibility is to fix the two 
parameters Kj and Kw modeling insolation and to identify only the remaining 
five which are then uniquely identifiable. Especially in the context of numerical 
tests, this is a convenient procedure since the model itself remains unaltered. Al¬ 
ternatively, by replacing the critical reaction terms altogether, new parameters 
without dependencies can be introduced. The analysis allowed us to formulate 
different substitutes in such a way that the reaction terms’ mathematical behav¬ 
ior remains similar to the original one. Other alternatives could be developed. 
As a future task, it remains to find out which alternative is in line with the 
ecosystem modeled. 

As a conclusion, the analysis improves the interpretation of tests in the 
context of parameter identification and the assessment of methods tested. By 
using one of our proposed alternatives instead of the original N-DOP model, 
one important source of uncertainty can be eliminated. 

Finally, the analysis concerning identifiability and solvability of the N-DOP- 
Fe model is performed with the help of universal mathematical methods. Thus, 
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the analysis conducted for the N-DOP-Fe model provides the basis for the 
assessment of other models as well. 
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Appendix A. Derivation of a characterization for dependencies 


Oi i - 


In order to find a characterization for dependencies, we solve the equation 

lit ii-t T p x%Kw 2 

y*- r r„ — xzKwi I xnKwi I ■ ts _ c 


■{/e"* 3Kwi |/e 


Kll = «21 


ls/i| + Kpi 1 ~ \vi\ + Kp 2 \Ie~*» K v”\ + K I2 
for |yi|. Taking into account ai,Ie~ XsKwi > 0 and y-\ ^ 0, we obtain 

a 1 _ |2/i| + Kpi le~ XsKwi + K n yi i e -x 3 K W2 

a 2 ~ yi i e -x 3 i<wi \ yi \ +K P2 Ie ~ X3Kw 2 +1< 12 

= |yi| + Kpi Ie- X ^ K ^ + K n c - x „ (Kw ,- Kw ,) 

\yi\ + Kp 2 Ie ~ xaKw2 + Kj2 

Using the abbreviation C := 01 / 02 , we calculate 


0 = C(|yi| + K P2 ){Ie - X3Kw2 + A /2 ) - (|t/i| + K P1 )(Ie~ X3Kwi + K n )e~ x ^ Kw2 - Kwi) 
= C\ yi \{Ie~ X3Kw2 + K n ) + CKp 2 {Ie~ X3Kw2 + A /2 ) 

- lyiK/e"® 3 ^ 2 + K I1 e- X3( ' Kw2 ~ Kwi) ) - K P1 (Ie~ X3Kw2 + K n e ~ X3(Kw2 - Kwi) ) 
= |yi|{C , (/e" a;3 * :iV2 + K I2 ) - {Ie ~ XaKw2 +A' n e“ l3(%2 - %l) )} 

+ CK P2 {Ie- X3Kw2 + K I2 ) - K P1 (Ie~ X3Kw2 + K n e- X3{Kw2 - Kwi) ). 


Rearranging the summands leads to 

\yi\{((C - l)I-K n e X3Kwi )e~ X3Kw2 + CI< I2 } 

= {(K P 1 - CKp 2 )I + K P1 Kne X3Kwi )e~ X3Kw2 - CK P2 K I2 . 

Taking into account Kj ^ 0, it is possible to prove that dependencies exist if 
and only if the expression in curly brackets is not zero. Dividing the equation 
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by this expression and canceling out e XsKw2 on the right-hand side, we obtain 

(K P i - CK P2 )I + K P1 K n e X3Kw ' - CK P2 K I2 e^ K ^ 

12/11 (C-1)1 - Kne^m+CK I2 e x 3 K w> ' [ ’’ 

Equation (IA.1I) corresponds to the condition specified in Sec. ED The additional 

condition for the constants ci,..., eg reflects the dependence of the coefficients 

in Eq. (IA.1I) on each other. 
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